clc; clear; close all;
cd 'data'

lambda_data = csvread('lambda.csv');
rho_data = csvread('rho.csv');
Y_data = csvread('Y.csv');

cd ..

I = 44;

alpha = 1; 
theta = 2.5; sigma = 3;
YEAR = 15;

for t = 1:YEAR
   lambda{t} = reshape(lambda_data(:,t),I,I);    
   rho{t} = 1-reshape(rho_data(:,t),I,I);    
   Y(:,t) = Y_data(:,t);
   rho_autarky(:,t) = max(0,1-(1-diag(rho{t}))./diag(lambda{t}));
end

for t = 1:YEAR
   LL(:,t) = diag(lambda{t});    
end

Real_change  = LL.^(-1/theta);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t = YEAR;

agrid = [1, 0.9, 0.8, 0.7 ,0.6, 0.5, 0.4 ,0.3, 0.2,0.1];

for a = 1:length(agrid)

alpha = agrid(a);    
    
s = 0.95; star = rho{t};

Lorenz5(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*sum((ones(I,I).*(s>star))...
    .*lambda{t}.*repmat(Y(:,t)',I,1)./repmat(Y(:,t),1,I).*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho{t})).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho{t})),2);

Share5(:,a) = 1 - Lorenz5(:,t);
 
s = 0.95; star_autarky = rho_autarky(:,t);
Lorenz5_autarky(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*(s>star_autarky)...
   .*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho_autarky(:,t))).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho_autarky(:,t)));
Share5_autarky(:,a) = 1 - Lorenz5_autarky(:,t);

end

Table = Share5.\Share5_autarky;

clearvars -except Table




